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We investigate the vortex state with Fulde-Ferrell-Larkin-Ovchinnikov (FFLO) modulations sug- 
gested for a high field phase of CeCoIns. On the basis of the quasiclassical Eilenberger theory, we 
calculate the three dimensional structure of pair potentials, internal magnetic fields, paramagnetic 
moments, and electronic states, for the s-wave and the d-wave pairings comparatively. The 7r-phase 
shift of the pair potential at the FFLO nodal plane or at the vortex core induces sharp peak states 
in the local density of states, and enhances the local paramagnetic moment. We also discuss the 
NMR spectrum and the neutron scattering as methods to detect the FFLO structure. 
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I. INTRODUCTION 

The Fulde-Ferrell-Larkin-Ovchinnikov (FFLO) stated 
is an exotic superconducting state expected to appear at 
low temperatures and high fields, when the paramagnetic 
depairing is significant due to the Zeeman effect under a 
magnetic field. In the FFLO state, since the Fermi sur- 
faces for up-spin and down-spin electron bands are split 
by the Zeeman effect, Cooper pairs of up- and down-spins 
acquire non-zero momentum for the center of mass coor- 
dinate of the Cooper pair, inducing the spatial modula- 
tion of the pair potential,M£&L&2, The possible FFLO 
state is widely discussed in various research fields, rang- 
ing from superconductors in condensed matter, neutral 
Fermion superfluids in an atomic cloud , 10 ' 11 i 12 i 13 to color 
superconductivity in high energy physics*^ 

Experimentally, the FFLO state is suggested in a 
high field phase of a quasi-two dimensional (Q2D) heavy 
Fermion superconductor CeCoInsi 15 : 16 Anomalous be- 
haviors of the sound wave velocity^ thermal conductiv- 
ity^ penetration depth^ and the NMR spectru m 20 : 21 
in the high field phase are considered as evidence of the 
FFLO structure. There, it is supposed that nodal planes 
of the pair potential run perpendicular to the vortex lines. 

In theoretical studies, many calculations for the FFLO 
state have been done by neglecting vortex structure. 
However, we have to include the vortex structure in ad- 
dition to the FFLO modulation, because the FFLO state 
appears at high fields in the mixed states.— Among the 
FFLO states, there are two possible spatial modulation 
of the pair potential A. One is the Fulde-Ferrell (FF) 
stated with phase modulation such as A oc e lC ^ z , where Q 
is the modulation vector of the FFLO states. The other 
is the Larkin-Ovchinnikov (LO) stated with the ampli- 
tude modulation such as A oc sinQz, where the the pair 
potential shows the periodic sign change, and A = at 
nodal planes. We discuss the case of the LO states in this 
paper, since some experimenta l 17 : 18 ' 20 and theoretical 7 ^ 
works support the LO state (at least in low temperature 
region) for the FFLO states in CeCoIn 5 . In the FFLO 
vortex state, it is instructive to analyze the role of the 



nodal plane in the LO states, which may give a clue to 
obtain clear evidence of the FFLO states among the ex- 
perimental data. 

When we consider vortex structure in the LO state, 
there are two possible choices of the configuration for 
the vortex lines and the FFLO modulation. That is, 
the modulation vector of the FFLO state is parallel^ or 
perpendicular— £ to the applied magnetic field. In our 
study, the former case is investigated by the quasiclassi- 
cal theory! ^ 22 : 23 : 24 : 25 In a previous quasiclassical study 
by Tachiki et aiA for this case, they calculate the struc- 
ture of the FFLO modulation after reducing the FFLO 
states with Abrikosov vortex lattice to a problem of the 
one dimensional (ID) system along the magnetic field di- 
rection. There, the radial dependence from the vortex 
core could not be analyzed, since the spatial dependence 
of vortex is approximately changed to the averaged one. 
In our study, fully three dimensional (3D) structures of 
the vortex and the FFLO modulation are determined by 
the selfconsistent calculation based on the quasiclassical 
theory. That is, we also selfconsistently obtain the radial 
dependence around vortex in the FFLO vortex states, 
which was not done before. 

On the other hand, the vortex and FFLO nodal plane 
structures in the FFLO state were calculated by the Bo- 
goliubov de-Gennes (BdG) theory for a single vortex in 
a superconductor under a cylindrical symmetry situa- 
tion^ This study clarifies that the topological structure 
of the pair potential plays important roles to determine 
the electronic structures in the FFLO vortex state. The 
pair potential has 27r-phase winding around the vortex 
line, and 7r-phase shift at the nodal plane of the FFLO 
modulation. These topologies of the pair potential struc- 
ture affect the distribution of paramagnetic moment and 
low energy electronic states inside the superconducting 
gap. For an example, the paramagnetic moment is en- 
hanced at the vortex core and the FFLO nodal plane. 
These structures are related to the bound states due to 
the 7r-phase shift of the pair potential. 

The other method to microscopically study the vortex 
state is the calculations by the quasiclassical Eilenberger 
theory. In the calculation by the BdG theory, we have to 
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assume a cylindrical symmetry and open boundary condi- 
tion for the system, because of the limit of matrix size to 
numerically solve the BdG equation in the present status 
of the computers. However, using the Eilenberger theory 
in the limit kpC ^ 1 (&F is Fermi wave number and £ is 
the superconducting coherence length), we are free from 
these restrictions. Therefore, we can discuss the cases 
of anisotropic pairing symmetry, such as <i-wave pairing, 
and arbitrary shape Fermi surfaces, where we obtain bro- 
ken cylindrical structures around vortices. Since we can 
consider the system of vortex lattice and periodic FFLO 
modulation by the periodic boundary condition, we can 
discuss the overlaps between tails of the neighbor vortex 
cores or FFLO nodal planes. These calculations for the 
periodic systems make us possible to estimate the res- 
onance line shapes in the NMR experiments and form 
factors in the neutron scattering experiment. 

The purpose of this paper is to investigate the 3D 
structures of the vortex and the FFLO nodal plane in 
the FFLO vortex state by the quasiclassical theory. We 
calculate the spatial structures of pair potentials, param- 
agnetic moments, internal magnetic fields and electronic 
states in the vortex lattice state under the given period of 
the FFLO modulation. From the results, we confirm the 
fundamental properties of the FFLO vortex states ob- 
tained by the BdG theory, and clarify the further details 
of the FFLO vortex states in the system of the vortex 
lattice and periodic FFLO modulation. Since the super- 
conducting pairing symmetry of CeCoIn 5 is suggested to 
be a e?-wave with line nodes, we calculate both cases of 
the isotropic s-wave pairing and the anisotropic d-wave 
paring comparatively, in order to study the contribution 
of the pairing symmetry to the FFLO vortex state^i 2 ! 
In the calculation, we use the Q2D Fermi surface and a 
magnetic field is applied to the ab plane, as is done for 
CeCoIn 5 . As possible experimental methods to detect 
the FFLO modulation, we discuss the NMR resonance 
line shaped and the neutron scattering experiments. We 
show that the spatial structure of the paramagnetic mo- 
ment reflects the NMR resonance line shape in the NMR 
experiment observing Knight shift under an applied mag- 
netic field. 

After introducing the formulation for the calculation of 
the FFLO vortex structure by the quasiclassical theory 
including the paramagnetic effect in Sec. [Ill we investi- 
gate the spatial structure of the FFLO vortex states in 
Sec. IIII1 and electronic states in Sec. IIVI Using the cal- 
culated structure of the FFLO vortex states, we discuss 
the NMR resonance line shape in Sec. IVland the neutron 
scattering in Sec. I VII The last section is devoted to the 
summary and discussions. 



II. QUASICLASSICAL THEORY INCLUDING 
PARAMAGNETIC EFFECT 

In this study for the FFLO vortex state, we take ac- 
count of the paramagnetic depairing effect due to the 



Zeeman splitting term [i B B(r) in addition to the orbital 
depairing effect by the vector potential A(r). Therefore, 
Hamiltonian for the spin-singlet pairing superconductor 
is given by 

U-^N= /d 3 r ^(r)^(r)^(r) 

-J d 3 n y^rajAfri.ra^fCri^Ifra) 

+A*(r 1 ,r 2 )V> i (r 2 )V T (ri)} (1) 

with 

K a {v) = ^(- + JaV + ^bSW-Mo, (2) 
2m \ l 0o J 

the flux quantum (f>o, and a = ±1 for up/down spin 
electrons. In the following, length, temperature, Fermi 
velocity, magnetic field and vector potential are, re- 
spectively, scaled by i?oi T c , vp, B and B R . Here, 
i?o = Tlv? /2-Kk B T c , B — frc/2\e\R 2 ,, and vf — (v 2 ,)]/ 2 is 
an averaged Fermi velocity on the Fermi surface. (■ • -)k 
indicates the Fermi surface average. Energy E, pair 
potential A and Matsubara frequency lui are scaled by 

Using the quasiclassical Green's functions g(uii + 
ip,B, k, r), /(tjj+i/LB, k, r), f(u>i+ijlB, k, r), Eilenberger 
equations are given by^ & 22 i 23 i 24 ' 25 

{lui + i£B + v • (V + iA)} / = A(r, k)g, 
{cjt + i/LB - v • (V - iA)} /t = A*(r, k)g, 
vV 5 = A*(r,k)/-A(r,k)/t, (3) 

where g = (1 - fpf 2 , Reg > 0, A(r,k) = A(r)0(k), 
and fi = fi B B /irk B T c . In Eq. ©, r = (r x + r 2 )/2 is 
the center of mass coordinate of the Cooper pair, and 
k = (k a , kb, k c ) is a relative momentum of the Cooper 
pair. We set the pairing function 0(k) = I in the s-wave 
pairing, and (f>(k) = x 2_ J/2 (k) = V2(fc 2 - fc 2 )/(fc 2 + k 2 ) 
or 4>(k) = (j) X y(k) = 2 v / 2fc a fcb/(fc 2 + k 2 ) in the ci-wave 
pairing. The vector potential is given byA= |Bxr + a 
in the symmetric gauge, with an average flux density B. 
The internal field is obtained as B(r) = B + V x a. 
The pair potential is selfconsistently calculated by 

A(r)=. 9o iV T V^ k )( / + /t *)) k W 

0<oJi<w cut 

with (goNo)- 1 = lnT + 2TE <^<^ ut ^" 1 - We use 
c<j cu t = 20k B T c . The vector potential is selfconsistently 
determined by the paramagnetic moment and the super- 
current as 

V x V x a(r) = j s (r) + V x M para (r) = j(r), (5) 
where 

j s (r) = -— ]T (v F Im{. 9 }) k (6) 
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is the supercurrent, and 

M para (r) = M ^1 - || ( Im {3»kj (7) 

is the paramagnetic moment in the vortex lattice state. 
Here, the normal state paramagnetic moment Mq = 
(p,/k) 2 B, k — B Q /TTkBT c ^/8nNo, N is the density of 
states at the Fermi energy in the normal state. 

When we calculate the electronic states, we solve Eq. 
([3|) with iuoi — > E+irj. The local density of states (LDOS) 
is given by 

}>k (8) 

for each spin component. We typically use r\ = 0.01. 

As a model of Fermi surface in CeCoIns, we use 
a Q2D Fermi surface with rippled cylinder-shape, and 
the Fermi velocity is given by vp = (v a > u 6> u c) °c 
(cos#, sin#, v z sinfc c ) at the Fermi surface kp = 
(k a ,kb,k c ) oc (fcpo cos#, fcpo sin#, fc c ) with — 7r < 9 < tt 
and — 7r < fc c < 7Ti^ In our calculation we set v z — 0.5, so 

that the anisotropy ratio 7 = £ c /£ Q b — ( w c)i/ 2 / (Oi/ 2 ~ 
0.5. A magnetic field is applied along the a axis direction 
in our calculation. Thus, the coordinate (x, y, z) for the 
vortex structure corresponds to (6, c, a) of the crystal co- 
ordinate. In the d-wave pairing, the case of the pairing 
function <f> x 2_ y 2 {4> xy ) corresponds to the configuration 
when B is applied to the anti-nodal direction (the nodal 
direction) of the superconducting gap. 

We solve Eq. © and Eqs. (4}- (7]) alternately, and ob- 
tain selfconsistent solutions as in previous works ; 23 ' 24 i 25 
by fixing a unit cell of the vortex lattice and a period 
L of the FFLO modulation. The unit cell of the vor- 
tex lattice is given by r = w\{vl\ — u 2 ) + w 2 u 2 + w 3 u 3 
with -0.5 < Wi < 0.5 (i=l, 2, 3), u x = (a, 0,0), u 2 = 
(a/2, ay, 0) and u 3 = (0,0, L). Reflecting anisotropy ra- 
tio 7, we set a y /a = v3t/2 with 7 = 0.5. For the FFLO 
modulation, we assume A(x,y,z) = A(x,y,z + L) and 
A(x, y, z) = -A(x, y, -z). Then, A(r) = at the FFLO 
nodal planes z = 0, and ±0.5i. These configurations 
of the FFLO vortex structure are schematically shown in 
Fig. [TJ which show the unit cell in the xz plane including 
vortex lines, and in the xy plane. 

We divide Wi to A^-mesh points in our numerical cal- 
culation, and obtain the quasiclassical Green's functions, 
A(r), M para (r) and j(r) at each mesh point in the 3D 
space. Typically we set N\ = N2 = N3 = 31. In the 
selfconsistent calculation of a, we solve Eq. ([SJ in the 
Fourier space q TOliTO2)m3 , taking account of the current 
conservation V • j(r) = 0, so that the average flux density 
per unit cell of the vortex lattice is kept constant. The 
wave number q is discretized as 

Vlmi ,m2 ,1113 

= TOiqi + W2qi + "i 3 q 3 (9) 

with integers m, (i — 1,2,3), where qi = 
(2ir/a,—ir/a y ,0), q 2 = (2Tr/a,ir/a y ,0), and 
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FIG. 1: Configurations of the vortex lines and the FFLO 
nodal planes are schematically presented in the xz plane in- 
cluding vortex lines (a) and in the xy plane (b). The vortex 
distance is a in the x direction, and the distance between the 
FFLO nodal planes is L/2. The hatched region indicates the 
unit cell. In (a) , along the trajectories presented by "0 — > tv" , 
the pair potential changes the sign (H > — ) across the vor- 
tex line or across the FFLO nodal plane, due to the 7r-phase 
shift of the pair potential. Along the trajectory presented by 
"0 — > 2n" , the sign of the the pair potential does not change 

(H > +) across the intersection point of the vortex line and 

the FFLO nodal plane, since the phase shift is 2n. In (b), • 
indicates the vortex center, ui — 112 and 112 are unit vectors 
of the vortex lattice. 

q 3 = (0, 0, 2-k/L). The lattice momentum is 
defined as G(q mi , m2iI „ 3 ) = (G x ,G y ,G z ) with 
G x = [JVisin(27rmi/iVi) + N 2 sm(2nm2/N 2 )]/a, 
G y = [-JVisin(27rmi/iVi) + N 2 sin{2TTm 2 /N 2 )}/2a y , 
and G z = N 3 sm(2iTm3/N 3 )/L. We obtain the Fourier 
component of a(r) as a(q) = j'(q)/|G| 2 , where 
j'(q) = j(q) — G(G-j(q)) /|G| 2 ensuring the current 
conservation V • j'(r) = 0, and j(q) is the Fourier 
component of j(r) in Eq. ([S])^ The final selfconsistent 
solution satisfies V • j(r) = 0. When we solve Eq. ^ by 
the explosion method, 2 — i 23 i 24 i 25 W e estimate A(r) and 
A(r) at arbitrary positions by the interpolation from 
their values at the mesh points, and by the periodic 
boundary condition of the unit cell including the phase 
factor due to the magnetic field. 



III. FFLO VORTEX STRUCTURE 

We consider the case of large GL parameter k = 
20 and large paramagnetic contribution Jl = 1.7. In 
these parameters, the upper critical field is suppressed 
to H c2 (T = 0.1T C ) - 0.22Bo, as a first order transi- 
tion. We discuss the FFLO vortex state when we set 
L = 50R or 70i? at T = 0.1T C and B = 0.2B , giv- 
ing Mo = 0.00139-Bo and the vortex lattice constant 
a = 8.5Rq. 

In the upper panels of Fig. [21 we show the spatial 
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(a) (b) (c) 

FIG. 2: (color online) Spatial structure of the FFLO vortex state in the xz plane at B — 0.2Bo, T = 0.1T C and L = 50-Ro- 
(a) Amplitude of the pair potential |A(r)|. (b) Paramagnetic moment M pa ra(r). (c) Internal magnetic field B z (r). The upper 
panels show the spatial variation within a unit cell [hatched region in Fig. QJa)] for the s-wave pairing. The middle panels 
present the profiles along the path UNCVU shown in Fig. [Ha) for the s-wave pairing. The lower panels are the same as the 
middle panels, but for the d-wave pairing. The solid (dashed) lines are for the case when the applied field is applied to the 
anti-nodal (nodal) direction. 



structure of the FFLO vortex state within a unit cell in 
the slice of the xz plane, i.e., the hatched region shown 
in Fig. QJa). In the middle and lower panels of Fig. 
the profiles of the spatial structure are presented along 
the path UNCVU shown in Fig. [Ha). The point C (x = 
y = z = 0) is the intersection point of a vortex and a 
nodal plane. The point V (x — y — and z = L/4) 
is at the vortex center and far from the FFLO nodal 
plane. The point N (x = a/2, y = z = 0) is at the 
FFLO nodal plane and outside of the vortex. The point 
U (x = a/2, y = 0, z = L/4) is far from both the vortex 
and the FFLO nodal plane. In Fig. [21 in addition to the 
s-wave pairing case in the middle panels, we also present 
the d-w&ve pairing cases in the lower panels in order to 
show that the qualitative features of the FFLO vortex 
structure are independent of the pairing symmetry and 
the field orientation within the ab plane. 

Figure [UJa) shows the amplitude of the order parame- 
ter. In the upper panel, we see that |A(r)| is suppressed 
near the vortex center at x — y = and the FFLO nodal 
plane at z — 0, ±0.5L. Far from the FFLO nodal plane 
such as z — 0.25L, |A(r)| shows a typical profile of the 
conventional vortex. When we cross the vortex line or 
the FFLO nodal plane, the sign of A(r) changes due to 
the 7r-phase shift of the pair potential as schematically 
shown in Fig. [IJa). In the profile of |A(r)| presented in 
the middle and the lower panels of Fig. Ufa), |A(r)| = 
along the FFLO nodal plane NC and along the vortex 
line CV. While amplitude |A(r)| in the s-wave pairing 
[middle panel] is larger than that in the d-wave pairing 
[lower panel], this is due to the difference of the ratio 



A(H = 0)/T c depending on the pairing symmetry. 

Correspondingly, paramagnetic moment M para (r)/Afo 
is presented in Fig. [2fb). As is well-known as the Knight 
shift, the paramagnetic moment is suppressed in the spin- 
singlet pairing superconducting state. In the figures, we 
see that M para (r) is suppressed outside of the vortex and 
far from the FFLO nodal plane, as expected. However, 
-^para(r) is enhanced at the vortex core or at the FFLO 
nodal plane, both in the s-wave and the d-wave pairings. 
The reason for these structures of M para (r) is discussed 
later in connection with the LDOS. At the FFLO nodal 
plane M para (r) ~ M [path NC in Fig. Hfb)]. Along the 
vortex line, M para (r) is enhanced more than Mq far from 
the FFLO nodal planes [position V in Fig. [2tb)]. 

Figure [2jc) presents the z-component of the internal 
field, B z (r). Due to the contribution of the enhanced 
-^para(r), B z (r) is enhanced at the FFLO nodal plane 
even outside of the vortex. A part of the contributions 
by M para (r) is compensated by the diamagnetic contri- 
bution, because the average flux density per unit cell of 
the vortex lattice in the xy plane should conserve along 
the magnetic field direction. Therefore, due to the con- 
servation, the enhancement of B z (r) at the FFLO nodal 
plane [path NC in Fig. 0Jc)] is smaller, compared with 
the enhancement of M para (r) at the FFLO nodal plane 
[path NC in Fig. Hfb)]. While B z (r) is largely enhanced 
than B at the vortex core far from the FFLO nodal plane 
[position V in Fig. [5Jc)], B z (r) is not largely enhanced 
at the vortex core in the FFLO nodal plane [position C] . 
Therefore B z (r) ~ B = 0.2 at the FFLO nodal plane 
[path NC]. 
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FIG. 3: (color online) Spectrum of the LDOS for up-spin elec- 
trons Nf (E, r) /No (solid lines) and for down-spin electrons 
Ni(E,t)/Nq (dashed lines) in the s-wave pairing. B — 0.2Bo, 
T — 0.1T C and L = 50-Ro- From top panel to bottom, we 
present the LDOS at positions U, V, N, and C. Their loca- 
tions are shown in Fig. []Ja). 
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FIG. 4: (color online) The same as Fig. [3] but in the d-wave 
pairing, (a) B is applied to the anti-nodal direction, (b) B is 
applied to the nodal direction. 



IV. ELECTRONIC STRUCTURE IN THE FFLO 
VORTEX STATE 



We compare the behaviors of the middle panels for the 
s-wave pairing and the lower panels for the d-wave pair- 
ing. The profile along the path VU is almost the same as 
that of the conventional vortex without FFLO modula- 
tion. Quantitatively, the variable ranges of M para (r) and 
B z (r) are smaller in the ci-wave pairing, compared with 
the s-wave pairing. In the lower panels, we also show the 
dependence on the relative angle of the magnetic field 
and the gap nodes. In any cases of pairings and the 
magnetic field directions, we see qualitatively the same 
behavior of the FFLO structure, i.e., M para (r) ~ Mq 
and B z (y) ~ B. Therefore, if we study the dependence 
on the pairing symmetry and the magnetic field direc- 
tion within the ab plane, we need to carefully examine 
the quantitative difference of the FFLO structure in the 
calculation of the pair potential, the paramagnetic mo- 
ment and the internal magnetic field. On the other hand, 
the dependences on the pairing symmetry appear more 
clearly in the structure of the electronic states in the vor- 
tex states^ Therefore, we discuss the electronic struc- 
ture in the FFLO vortex states in the next section. 



The LDOS spectrum for up- and down-spin electrons 
are presented at some positions in Fig. [3] for the s-wave 
pairing and in Fig. [¥]for the <i-wave pairing. In the qua- 
siclassical theory, N a (E,r) are symmetric by E <-> —E 
in the absence of the paramagnetic effect [Jx = 0). In 
the presence of the paramagnetic effect, the LDOS spec- 
trum for up- (down-) spin electrons is shifted to positive 
(negative) energy by flB due to the Zeeman shift. In this 
case, we have a relation N-[ (E, r) = Ni(—E, r) within the 
quasiclassical theory. 

Far from the FFLO nodal plane and outside of vortex, 
as shown in the spectrum of the top panel for the po- 
sition U in Fig. [31 we see the superconducting full-gap 
structure of the s-wave pairing. There, small LDOS also 
appears at low energies inside the gap due to the low en- 
ergy excitations extending from the vortex cores and the 
FFLO nodal planes at finite magnetic fields. Since the 
LDOS at E < (E > 0) are occupied (empty) states, 
there is a relation between the LDOS spectrum and local 
paramagnetic moment as 

Mpara(r) = -/x B / (Af-f (E, r) - iV 1 ( J B,r))d£. (10) 
In the full-gap case of the s-wave pairing, the difference of 
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the occupation number between up- and down-spin elec- 
trons is small, since the LDOS below the gap are occupied 
similarly in N^(E,r) and N±(E,r). This is the reason 
why M para (r) is suppressed at the position U. Small but 
finite M para (r) comes from the small LDOS weight of low 
energy states inside the gap in the top panel of Fig. [3] 

The LDOS spectra at the position V on the vortex 
center and at the position N on the FFLO nodal plane 
are, respectively, presented in the second and the third 
panels of Fig. [3] In these figures, N^(E, r) and Ni(E, r), 
respectively, have a sharp peak at E = n+ and E = 
with fi-j. = fiQ ± jlB. These peaks are related to the 
topological structure of the pair potential, as schemati- 
cally shown in Fig. [1] Since a vortex has phase winding 
2ir, along the trajectory through the vortex center, A(r) 
changes the sign by the 7r-phase shift across the vortex 
center. Also at the trajectory through the FFLO nodal 
plane, A(r) changes the sign across the nodal plane. The 
bound states appear when the pair potential has the ir- 
phase shift, and form the "zero-energy peak" . This peak 
is shifted to E = /i+ or E = /i_ due to the Zeeman 
effect i 26 i 28 Since the peak of the LDOS spectrum for up- 
spin electrons is an empty state (E > 0) and the peak 
of the LDOS for down-spin electrons is an occupied state 
(E < 0), M para (r) becomes large at these positions, from 
the relation in Eq. (fT0|) . 

On the other hand, along the trajectory through the 
intersection point of a vortex and a nodal plane, A(r) 
does not change the sign, because the phase shift is 2tt 
by summing ir due to vortex and ir due to the nodal plane, 
as schematically shown in Fig. [T] Thus, the sharp peaks 
do not appear at E = fi± as seen from the bottom panel 
of Fig. [3j Instead, JV-j- (E, r) has two broad peaks at finite 
energies shifted upper or lower from In this situation, 
-Mpara(r) is still large at position C, as in positions V and 
N, since the LDOS in both peaks are empty (E > 0) in 
N^(E,r), and occupied (E < 0) in N±(E,r). 

The differences of the s-wave and the <i-wave pairings 
appear in the LDOS spectrum of the superconducting 
gap. At a zero field, the LDOS has the full-gap structure 
in the s-wave pairing, while it has a V-shape gap in the d- 
wave pairing with line nodes. Under the magnetic field, 
these features of the gap structure are smeared by low 
energy excitations under the magnetic field, but can be 
seen far from the vortex and the FFLO nodal plane, as is 
shown in top panels [position U] in Figs. [3]and|U In the 
d-wave pairing, the LDOS spectrum has larger weight of 
low energy states within the gap [top panels in Fig. 0], 
and M para (r) is larger [position U in Fig. E{b)], compared 
with the s-wave pairing. From Fig. [51 we see that in the 
e?-wave pairing the LDOS structures at the vortex line 
and at the FFLO nodal plane are qualitatively the same 
as in the s-wave pairing. The sharp peak structure at 
E = n± appears in the LDOS at the vortex line [second 
panels] and at the FFLO nodal plane [third panels] . And 
the peak is split to two broad peaks, and shifted to upper 
or lower from fi± at the intersection point of the vortex 
line and the FFLO nodal plane [bottom panels for the 




(b) (c) 

FIG. 5: (color online) FFLO vortex structure at B — 0.15-Bo, 
T — 0.1T C and L = 50iio in the s-wave pairing, (a) Spectrum 
of the LDOS for up-spin electrons Nf(E,r)/No (solid lines) 
and for down-spin electrons N[ (E, r) /No (dashed lines) at the 
position C. (b) Paramagnetic moment M para (r). (c) Internal 
field -B z (r). The upper panels of (b) and (c) show the spatial 
variation within a unit cell. The lower panels present the 
profile along the path UNCVU. 



position C] . Figure 0] presents two cases of the magnetic 
field orientation, which show that the characteristic peak 
structures of the LDOS in the FFLO vortex structure 
do not significantly depend on the pairing symmetry and 
the relative angle of the magnetic field orientation and 
the gap-node direction. 

In order to understand the relation between the LDOS 
and the paramagnetic moment, it is instructive to see the 
lower magnetic field case at B = 0.15£>o shown in Fig. 
El Here, a = 9.8i? and M = 0.00105B . Figure [^a) 
shows the LDOS at the position C. Compared with the 
LDOS in Fig. [3j the Zeeman shift fiB of the spectrum is 
smaller. In this case, one of two broad peaks is occupied 
(E < 0) and the other is empty (E > 0) both for N^(E, r) 
and N±(E,r). Therefore, M para (r) is suppressed at the 
position C, i.e., the intersection point of the vortex line 
and the FFLO nodal plane, as shown in Fig. Ob), while 
-Mpara(r) is larger at the vortex line (position V) and 
at the FFLO nodal plane (position N). This structure 
corresponds to that obtained in the calculation by the 
BdG theory^ Also in the spatial structure of the internal 
field, B z (r) is suppressed at the position C as shown in 
Fig. [5)^c) by the contribution from the spatial structure 
of Af para (r). 

In Fig. [6l we show the spectral evolution of 
Nf(E,r)/No around vortex and the FFLO nodal plane, 
when the pairing is <i-wave and magnetic field is applied 
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FIG. 6: (color) Spectral evolution of the LDOS for up- 
spin electrons along the path CVUNC shown in Fig. Ola). 
Density plot (a) and stereographic view (b) of iVj (E, r)/No 
are presented in the d-wave pairing, when magnetic field 
B = 0.2£?o is applied to the anti-nodal direction. T = 0.1T C 
and L — 50-Ro- 




FIG. 7: (color online) (a) Distribution function of the para- 
magnetic moment. We show P(M) as a function of M para . 
(b) Distribution function of the internal magnetic field. We 
show P(B) as a function of B z . The solid (dashed) lines are 
for the vortex state with (without) the FFLO modulation. 
B = 0.2B , T = O.lTc and L = 50R . The top panels are 
for the s-wave pairing. The middle (bottom) panels are for 
the d-wave pairing when the magnetic field is applied to the 
anti-nodal (nodal) directions. The heights of P(M) and P(B) 
are scaled so that J P(M)dM = J P(B)dB = 1. 



to the anti-nodal direction within ab plain. We obtain 
similar spectral evolution also in the s-wave pairing. This 
spectral evolution is compared with that obtaine d by 
BdG theory in the s-wave pairing [see Fig. 3 in Ref. l2q . 
The spectral structure is qualitatively good accordance 
with that. In the path VU, we see a typical spectral evo- 
lution of Caroli-de Gennes-Matricon states around the 
conventional vortex corei 28 ' 29 i 30 ' 31 i 32 With aparting from 
the vortex center, vortex core states have higher angular 
momentum, and the energy levels are shifted to higher 
energy. Along the vortex line CV, moving from the nodal 
plane, two peak states at C are merged to the conven- 
tional vortex core states with E = //+ at V. Also along 
the FFLO nodal plane NC, two peak states at the vortex 
center C are merged to the bound states with E = fi + at 
the midpoint N of two vortices. The calculation for the 
midpoint was not done in the previous BdG calculation 
assuming cylindrical symmetry^ 

V. NMR SPECTRUM 

In the NMR experiment, resonance frequency spec- 
trum of the nuclear spin resonance is determined by the 
internal magnetic field and the hyperfine coupling to the 
spin of the conduction electrons. Therefore, in a simple 
consideration, the effective field for the nuclear spin is 



given by B cS (r) = B(r) + A hi M paTa (r), where A M is a 
hyperfine coupling constant depending on species of the 
nuclear spins. The resonance line shape of NMR is given 
by 

P(u) = J 5{w - £ eff (r))dr, (11) 

i.e., the intensity at each resonance frequency u) comes 
from the volume satisfying uj = B e g(r) in a unit cell. 
When the contribution of the hyperfine coupling is dom- 
inant, the NMR signal selectively detects M para (r). This 
is the experiment observing the Knight shift in super- 
conductors^ 3 ^ As the NMR spectrum of the Knight 
shift, we calculate the distribution function P(M) = 
J S(M — M para (r))dr from the spatial structure of the 
paramagnetic moment M para (r) shown in Fig. [^b). We 
show P(M) in Fig. [7{a). On the other hand, in the 
case of negligible hyperfine coupling, the NMR signal is 
determined by the internal magnetic field distribution. 
This resonance line shape is called "Redfield pattern" of 
the vortex lattice. In Fig. [7£b), we show the distribu- 
tion function P(B) = J 5(B — _B z (r))dr calculated from 
the internal field B z (r). In Fig. \7\ the s-wave pairing 
case and two cases of the field orientations in the d-wave 
pairing are displayed. 

First we discuss the line shape of the distribution func- 
tion P(M). The spectrum of P(M) in the conventional 
vortex state without the FFLO modulation is shown by 



dashed lines in Fig. \7^a). There, the peak of P(M) comes 
from the signal from the outside of the vortex core. The 
temperature dependence of this peak position M gives 
Knight shift observed in the NMR experiments! 33 ! 34 In 
our calculation, we confirm that M is decreased from 
Mq on lowering temperature below T c . The spectrum 
of P(M) has a tail toward larger M by the vortex core 
contribution with large M para (r). The vortex core con- 
tribution is a ID structure, their volume contribution is 
small in the spectrum, compared with the peak inten- 
sity due to the large volume contribution by the outside 
of the vortex. While P{M) is slightly enhanced at the 
upper edge of M, this is the feature of the large param- 
agnetic effect. At the higher field near H C 2, the vortex 
core radius is expanded by the contribution of M para (r) 
enhanced at the vortex core, and the profile of M para (r) 
becomes flatter around the expanded vortex core region. 
These contributions slightly enhance P(M) at the upper 
edge. This enhancement is not seen when the paramag- 
netic effect is not large, such as at lower fields. 

In the vortex state with the FFLO modulation, the line 
shape P(M) becomes double peak structure, as presented 
by solid lines in Fig. [TJa). The height of the main peak 
decreases, and there appears a new peak coming from the 
FFLO nodal plane near M para ~ Mq. This is because a 
part of contribution from the outside region of the vortex 
is shifted to the contribution of the FFLO nodal plane. 
The contribution from the two dimensional structure of 
the FFLO nodal plane appears in P(M) more clearly 
than that of the ID structure of the vortex line. 

Second, we discuss the distribution function P(B) of 
the internal magnetic field presented in Fig. \T[b). There, 
in the absence of the FFLO modulation (dashed lines), 
the Rcdfield pattern P(B) has sharp peak corresponding 
the saddle point of the internal field distribution. The 
tail to higher B comes from the increasing magnetic field 
around the vortex core. The small enhancement of P{B) 
at the upper edge of B is also due to the contribution 
from the flat profile of M para (r) around the expanded 
vortex core, mentioned above. In the presence of the 
FFLO modulation (solid lines), the height of the original 
peak is decreased, and a new peak appears at B ~ B = 
0.2i?o as the contribution of the FFLO nodal plane. As 
a feature of P{B), the new peak appears near the main 
peak, compared with the line shape of P(M). 

In the ci-wave pairing seen in lower two panels in Fig. 
[7J while the ranges for the variation of M and B are 
narrower, the double peak structures of the FFLO vor- 
tex state are seen as in the s-wave pairing. We also see 
that the shapes of P(M) and P(B) are not significantly 
affected by the orientation of the applied magnetic field 
within the ab plane in the line node case of the <i-wave 
pairing. 

Next, we discuss how the spectra of P(M) and P(B) 
depend on the period L of the FFLO modulation, and 
on the temperature T. The L-dependences of P(M) and 
P(B) are shown in the upper panels of Fig. [HJ with the 
profiles of M para (r) and B z (r) in the lower panels. With 



P(M) 



P(B) 







L=70 













M 




















''30 \ 


u "' 




u u 







V v 






if \ L470 






/ ;30 


u u 




ti" IT 



(a) 



(b) 



FIG. 8: (color online) Structure of paramagnetic moment (a) 
and internal magnetic field (b) when L/Ro = 70 (solid lines) 
and 30 (dashed lines) in the d-wave pairing when the magnetic 
field is applied to the anti-nodal direction. B = 0.2£?o and 
T/T c = 0.1. The upper panels show the distribution function 
P(M) or P(B). The lower panels show the profile of M pa ra(r) 
or B z (r) along the path UNCVU. 
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FIG. 9: (color online) Structure of paramagnetic moment (a) 
and internal magnetic field (b) at T/T c =0.1 (solid lines) and 
0.2 (dashed lines) in the d-wave pairing when the magnetic 
field is applied to the anti-nodal direction. B = 0.2Bo and 
L = 70-Ro- The upper panels show the distribution function 
P{M) or P{B). The lower panels show the profile of M para (r) 
or B z (r) along the path UNCVU. 



decreasing L, the new peak of the FFLO nodal plane 
is enhanced both in P(M) and P(B), since the relative 
weight of volume contribution by the FFLO nodal plane 
increases in the distribution functions. As is seen from 
the profile along the path UN in the lower panel of Fig. [5J 
the effective width of FFLO nodal plane where M para (r) 
is enhanced near the FFLO nodal plane [position N] does 
not significantly change. Instead, the width of the region 
[around position U] with low M para (r) and B z (r), out- 
side of the effective FFLO nodal plane width, becomes 
shorter, when L decreases. 

The temperature dependences of P(M) and P{B) are 
shown in the upper panels of Fig. [9l with the profiles of 
-Mpara(r) and B z (r) in the lower panels. In the spectra 
of P(M) and P(B), the height of the new peak due to 
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the FFLO nodal plane becomes higher with increasing T. 
This is because the effective width of FFLO nodal plane 
becomes larger at higher T, since slope of M para (r) and 
B z (r) from the FFLO nodal plane decreases as shown in 
the lower panel of Fig. [9] [see slope around the position N 
in the path UN] . This enhancement of the new peak of the 
FFLO nodal plane at higher T comes from the increase 
of the effective coherence length. On the other hand, 
the slight enhancement of P(M) and P(B) at the upper 
edge, which appears even in the absence of the FFLO 
modulation, is smeared at higher T. This difference of 
the T-dependence can be used to distinguish the origin 
of the new peak observed in experiments. If the new 
peak becomes eminent at higher T, the origin of the new 
peak is the FFLO nodal plane. If the new peak becomes 
smeared at higher T, it may be an enhancement of the 
upper edge due to the paramagnetic vortex core effect. 
We also see that the main peak shifts to higher M in 
P(M) with increasing T, since M para (r) outside of vortex 
[around position U] is shifted to higher. 

The experimental observation of the NMR resonance 
line shape was used to identify the FFLO phase in 
CeCoIns^ There, the NMR spectrum shows the double 
peak structure in the FFLO phase, appearing new peak 
in addition to the main peak in the vortex state. This 
new peak is considered as a contribution of the FFLO 
nodal plane. This double peak structure is qualitatively 
reproduced by our calculation. In this status, our results 
may be qualitatively compared to the experimental data, 
since we did not perform fine tunings of parameters such 
as L, for the reason that we need further heavy calcula- 
tions. The mechanism of rf-field penetration in the FFLO 
stated may further modify the NMR spectrum from the 
original P{M). In the experimental data, the height of 
the new peak increases at higher T in the NMR spec- 
trum. 20 This feature is qualitatively consistent with the 
T-dependence in our results shown in Fig. O while it is an 
open question how L increases as a function of T. It in- 
dicates that, among the T- and L-dependences discussed 
in Figs. [8] and [9l the contribution of the T-dependence 
seems to be dominant than the //-dependence in the in- 
creasing T process. 



VI. NEUTRON SCATTERING 

The modulation of the internal magnetic field B z (r) 
may be observed by the neutron scattering. If the pe- 
riodic modulation along the z-direction is observed, it 
can be direct evidence of the FFLO modulation. There- 
fore, we discuss the neutron scattering in the FFLO 
vortex state. The intensity of the (h, k, ^-diffraction 
peak is given by I h ,k,i = \Fh,k,i\ 2 /\lh,k,i\ with the 
wave vector Qh,k,l given in Eq. Here we write 

(mi, 7712,771,3) = (h,k,l) following notations of the neu- 
tron scattering. The Fourier component Fh,k,i is given 
by B z (r) = J2h,k,l Fh > k > 1 ex P( ic Lh,M ' r )- In Fi g- EHa), we 
present /-dependence of \Fh,k,i\ 2 f° r Qi, k) — (1, 0), (1, 1) 
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FIG. 10: Factor |-F\,fc,i| 2 for the intensity of the neutron 
scattering, calculated from the internal magnetic field B z (r) 
(a) or from the paramagnetic moment M para (r) (b) in the 
d-wave pairing when the magnetic field is applied to the anti- 
nodal direction. |fh,fc,i| 2 is plotted as a function of I for 
(h,k) = (1,0), (1,1) and (2,0). In (b), we also show \F h , k ,l\ 2 
for (h,k) = (0,0). B = 0.2B , T = 0.1T C and L = 70R . 
Since the period of the FFLO modulation is L/2 for B z (r) 
and M para (r) along the z direction, the diffraction peaks ap- 
pear at spots with even I. 



and (2,0). The spot at (h,k,l) = (1,0,0) and (0,1,0) is 
used to determine the configuration and the orientation 
of the vortex lattice in the experiment of the small angle 
neutron scattering (SANS), and the higher component 
Fh.k.o is used to estimate the structure of the internal 
magnetic field B z (r)££- It is noted that Fo,o,o = B and 
Fo,o,l — for I ^ 0, because average flux density B within 
the unit cell of the vortex lattice is constant along the z- 
direction. Therefore, to detect the FFLO modulation, we 
have to use the spot (1, 0, 2), which has largest intensity 
among the spots related to the FFLO modulation. The 
spot (1, 0, 2) is next to the spot (1, 0, 0), which is used in 
the conventional SANS experiment to observe the stable 
vortex lattice configuration. 

If the spins of the paramagnetic moment M para (r) 
strongly couple to the lattice modulation, it is possi- 
ble to observe the structure of the paramagnetic moment 
through the lattice modulation by the neutron scattering. 
In this case, the Fourier component of the neutron scat- 
tering is given by M para (r) = J2h,k,i F h,k,i exp(iq^ jfe! ; • r). 
In Fig. [TUf bh we present \Fh^,i\ 2 calculated from 
•^para(r)- Since Fo,o,z (I 7^ 0) is not zero in this case, 
the intensity at (0,0, 1) is larger than that at (1,0, 1) for 
the same /. 



VII. SUMMARY AND DISCUSSION 

We investigated roles of the vortex and the FFLO 
nodal plane in the FFLO state, based on the quasiclas- 
sical Eilenberger theory. We selfconsistently calculated 
fully 3D spatial structure of the pair potential, the inter- 
nal magnetic field, the paramagnetic moment, and local 
electronic states in the vortex lattice state under given 
period of the FFLO modulation L. We have seen that the 
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topological structure of the pair potential determines the 
qualitative structure in the FFLO vortex state. At the 
FFLO nodal plane or at the vortex core, the 7r-phase shift 
of the pair potential gives rise to the sharp peaks in the 
local density of states, and enhances the local paramag- 
netic moment. Based on these results, we also discuss the 
experiments of NMR spectrum and neutron scattering, to 
identify characteristic behaviors in the FFLO states. 

In this work, we calculate the FFLO vortex structure 
both in the isotropic s-wave paring and in the anisotropic 
d-wave pairing with line nodes, when the FFLO nodal 
planes are perpendicular to the vortex direction. If we 
compare the FFLO vortex states under the same FFLO 
modulation period, we see that even in the d-wave pair- 
ing the qualitative structures of the FFLO vortex state 
are the same as in the s-wave pairing, and do not signif- 
icantly depend on the relative angle between the applied 
magnetic field direction and the node direction within 
the ab plane. This is because qualitative features re- 
ported in this paper come from the topological structure 
of A(r). Quantitatively, the dependences of the FFLO 
vortex state on the pairing symmetry or the magnetic 
field direction relative to the node of the superconduct- 
ing gap are important features to characterize the FFLO 
vortex state d 36 ' 37 For the purpose, we need estimates of 
the optimized FFLO period L and the phase diagram of 
the FFLO state in the mixed state. These calculations 
belong to future works within our framework of the self- 
consistent calculations for the 3D structure of the FFLO 
vortex states. However, qualitative features do not signif- 
icantly depend on this process. We did some calculations 
for different L, and confirmed that our results reported in 
this paper are qualitatively unchanged as for the FFLO 
vortex structure, because they come from the topological 
structure of the FFLO states, as discussed above. 

Compared to results in the previous work by BdG the- 
ory^ we obtain qualitatively the same features of the 
paramagnetic moment and local electronic states around 
vortices and FFLO nodal planes in the quasiclassical cal- 
culation. We note that these qualitative features by the 
7r-phase shift of the pair potential were confirmed in this 
study also for the cases of d-wave pairing and the Q2D 



Fermi surface with parallel applied fields along ab plane, 
which is not considered in the previous work by BdG the- 
ory assuming cylindrical symmetry structure and open 
boundary conditions. In this quasiclassical calculation 
with periodic boundary condition for the vortex lattice 
and FFLO modulation, we can take account of contri- 
butions by overlapping of electronic states from neighbor 
vortices or FFLO nodal planes. For example, we can 
discuss local structure at the mid-point between neigh- 
bor vortices. Compared with the previous quasiclassical 
work,— since we calculate fully 3D structure of the FFLO 
vortex states without mapping to ID problem, we can 
explicitly analyze the contribution of the vortex lines in 
the FFLO vortex states. Thus, we can quantitatively 
estimate the spatial structure around vortices along the 
radial direction, including the vortex structure within the 
FFLO nodal plane. There, we see exotic structure at the 
intersection point of vortex and FFLO nodal plane. By 
this 3D calculation of the FFLO vortex states, we can 
evaluate the NMR spectrum and form factors of neutron 
scattering, including vortex contributions in the FFLO 
states. 

From the spatial structure of the paramagnetic mo- 
ment in the FFLO vortex state, we can estimate the 
NMR spectrum in the Knight shift experiment. In the 
FFLO state, we obtain double peak structure of the NMR 
spectrum, as observed in the FFLO state of CeCoIns^ 
and discuss the dependences on the temperature and on 
the FFLO period. From the spatial structure of the in- 
ternal magnetic field, we also discuss the neutron scatter- 
ing in the vortex states as methods to detect the FFLO 
structure. Since this is a direct observation of the FFLO 
modulation, we hope that this feature will be observed 
in experiments to confirm the FFLO states. 
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